﻿#include "geom_surface_bspline.h"

CAGD::GeomBSplineSurface::GeomBSplineSurface(int udeg, int vdeg, const std::vector<std::vector<BPnt>> &cpts)
    : m_cpts(cpts)
{
    int p = udeg;
    int q = vdeg;
    int n = cpts.size() - 1;
    int m = cpts.front().size() - 1;

    // 生成均匀节点向量
    for (int i = 0; i < p + 1; i++)
        m_uknots.push_back(0);
    for (int i = p + 1; i < n + 1; i++)
        m_uknots.push_back(1.0 * (i - p) / (n + 1 - p));
    for (int i = n + 1; i < n + p + 2; i++)
        m_uknots.push_back(1);

    for (int i = 0; i < q + 1; i++)
        m_vknots.push_back(0);
    for (int i = q + 1; i < m + 1; i++)
        m_vknots.push_back(1.0 * (i - q) / (m + 1 - q));
    for (int i = m + 1; i < m + q + 2; i++)
        m_vknots.push_back(1);
}

CAGD::GeomBSplineSurface::GeomBSplineSurface(const std::vector<double> &uknots, const std::vector<double> &vknots,
                                             const std::vector<std::vector<BPnt>> &cpts)
    : m_uknots(uknots), m_vknots(vknots), m_cpts(cpts)
{
}

int CAGD::GeomBSplineSurface::NbUKnots() const
{
    return m_uknots.size();
}

const std::vector<double> &CAGD::GeomBSplineSurface::UKnots() const
{
    return m_uknots;
}

void CAGD::GeomBSplineSurface::SetUKnots(const std::vector<double> &knots)
{
    m_uknots = knots;
}

int CAGD::GeomBSplineSurface::NbVKnots() const
{
    return m_vknots.size();
}

const std::vector<double> &CAGD::GeomBSplineSurface::VKnots() const
{
    return m_vknots;
}

void CAGD::GeomBSplineSurface::SetVKnots(const std::vector<double> &knots)
{
    m_vknots = knots;
}

const std::vector<std::vector<CAGD::BPnt>> &CAGD::GeomBSplineSurface::Poles() const
{
    return m_cpts;
}

void CAGD::GeomBSplineSurface::SetPoles(const std::vector<std::vector<BPnt>> &cpts)
{
    m_cpts = cpts;
}

int CAGD::GeomBSplineSurface::UDegree() const
{
    return m_uknots.size() - m_cpts.size() - 1;
}

int CAGD::GeomBSplineSurface::VDegree() const
{
    return m_vknots.size() - m_cpts.front().size() - 1;
}

CAGD::GeomBSplineSurface::pnt_array_4d CAGD::GeomBSplineSurface::DerivCpts(int d, int r1, int r2, int s1, int s2) const
{
    int p = UDegree();
    int q = VDegree();
    int du = std::min(d, p);
    int dv = std::min(d, q);
    int r = r2 - r1;
    int s = s2 - s1;

    // 保存 P_{i,j}^(k,l) 对应 PKL[k][l][i][j]
    pnt_array_4d PKL(du + 1, pnt_array_3d(dv + 1, pnt_array_2d(r + 1, pnt_array(s + 1))));

    // 先处理 u 方向，得到 P_{i,j}^(k,0) 即 PKL[k][0][i][j]
    // 对 j 即 v 方向分别计算 u 方向的 k 阶导控制点
    for (int j = 0; j < s + 1; j++)
    {
        pnt_array_2d PK(du + 1, pnt_array(r + 1));

        // 0 阶导控制点（初始控制点） P_i^0
        for (int i = 0; i < r + 1; i++)
            PK[0][i] = m_cpts[r1 + i][s1 + j];

        // 计算 P_i^k
        for (int k = 1; k < du + 1; k++)
        {
            int tmp = p - k + 1;
            // P_i^k = (P_{i+1}^{k-1}-P_i^{k-1}) / (u_{i+p+1}-u_{i+k}) * (p-k+1)
            for (int i = 0; i < r - k + 1; i++)
                PK[k][i] = (PK[k - 1][i + 1] - PK[k - 1][i]) / (m_uknots[r1 + i + p + 1] - m_uknots[r1 + i + k]) * tmp;
        }

        // 储存计算结果，这部分只和 k,i 有关
        for (int k = 0; k < du + 1; k++)
            for (int i = 0; i < r - k + 1; i++)
                PKL[k][0][i][j] = PK[k][i];
    }

    // 再处理 v 方向，对于不同的 P_{i,j}^(k,0) 中的 k 计算 v 方向导
    for (int k = 0; k < du; k++)
    {
        // 对不同的 k，在 u 方向上计算出的控制点数为 r-k+1，因此对 i 从 0 到 r-k 操作
        for (int i = 0; i < r - k + 1; i++)
        {
            // u,v 方向导 不超过 d，这里减去已经对 u 求的 k 次导
            int dd = std::min(d - k, dv);
            pnt_array_2d PK(dd + 1, pnt_array(s + 1));

            // 0 阶导控制点（初始控制点），即从 PKL 中获得 P_{i,j}^(k,0) 控制点
            for (int j = 0; j < s + 1; j++)
                PK[0][j] = PKL[k][0][i][j];

            // 计算 P_{i,j}^(k,l)
            for (int l = 1; l < dd + 1; l++)
            {
                int tmp = q - l + 1;
                for (int j = 0; j < s - l + 1; j++)
                    PK[l][j] =
                        (PK[l - 1][j + 1] - PK[l - 1][j]) / (m_vknots[s1 + j + q + 1] - m_vknots[s1 + j + l]) * tmp;
            }

            // 最后获得控制点 P_{i,j}^(k,l)，这里只和 l,j 有关
            for (int l = 0; l < dd + 1; l++)
                for (int j = 0; j < s - l + 1; j++)
                    PKL[k][l][i][j] = PK[l][j];
        }
    }
    return PKL;
}

CAGD::GeomBSplineSurface::pnt_array_2d CAGD::GeomBSplineSurface::DerivPoints(int d, double u, double v) const
{
    pnt_array_2d SKL(d + 1, pnt_array(d + 1));

    int p = UDegree();
    int q = VDegree();
    int uSpan = Spline::FindSpan(u, p, m_uknots);
    int vSpan = Spline::FindSpan(v, q, m_vknots);
    int du = std::min(d, p);
    int dv = std::min(d, q);

    auto Nu = Spline::DerivAllBasis(uSpan, u, p, m_uknots);
    auto Nv = Spline::DerivAllBasis(vSpan, v, q, m_vknots);
    auto PKL = DerivCpts(d, uSpan - p, uSpan, vSpan - q, vSpan);

    for (int k = 0; k < du + 1; k++)
    {
        int dd = std::min(d - k, dv);
        for (int l = 0; l < dd + 1; l++)
        {
            for (int i = 0; i < q - l + 1; i++)
            {
                // N[j][r] = N_{r,j}
                BVec tmp;
                for (int j = 0; j < p - k + 1; j++)
                    tmp = tmp + PKL[k][l][j][i] * Nu[p - k][j];
                SKL[k][l] = SKL[k][l] + tmp * Nv[q - l][i];
            }
        }
    }
    return SKL;
}

void CAGD::GeomBSplineSurface::Bounds(double &U1, double &U2, double &V1, double &V2) const
{
    U1 = m_uknots[UDegree()];
    U2 = m_uknots[m_uknots.size() - UDegree() - 1];
    V1 = m_vknots[VDegree()];
    V2 = m_vknots[m_vknots.size() - VDegree() - 1];
}

void CAGD::GeomBSplineSurface::D0(double u, double v, BPnt &P) const
{
    int p = UDegree();
    int q = VDegree();
    int uspan = Spline::FindSpan(u, p, m_uknots);
    int vspan = Spline::FindSpan(v, q, m_vknots);
    auto n = Spline::Basis(uspan, u, p, m_uknots);
    auto m = Spline::Basis(vspan, v, q, m_vknots);

    BVec S;
    for (int i = 0; i < q + 1; i++)
    {
        BVec tmp;
        for (int j = 0; j < p + 1; j++)
            tmp = tmp + m_cpts[uspan - p + j][vspan - q + i] * n[j];
        S = S + tmp * m[i];
    }
    P = S;
}

void CAGD::GeomBSplineSurface::D1(double u, double v, BPnt &P, BVec &DU, BVec &DV) const
{
    auto DS = DerivPoints(1, u, v);
    P = DS[0][0];
    DU = DS[1][0];
    DV = DS[0][1];
}

void CAGD::GeomBSplineSurface::D2(double u, double v, BPnt &P, BVec &DU, BVec &DV, BVec &DUU, BVec &DVV,
                                  BVec &DUV) const
{
    auto DS = DerivPoints(2, u, v);
    P = DS[0][0];
    DU = DS[1][0];
    DV = DS[0][1];
    DUU = DS[2][0];
    DVV = DS[0][2];
    DUV = DS[1][1];
}

void CAGD::GeomBSplineSurface::DN(int N, double u, double v, std::vector<std::vector<BPnt>> &V) const
{
    V = DerivPoints(N, u, v);
}

std::vector<double> CAGD::GeomBSplineSurface::IsoUKnots() const
{
    int s = UDegree(), e = m_uknots.size() - UDegree();
    std::vector<double> knots;
    for (int i = s; i < e; i++)
        knots.push_back(m_uknots[i]);
    return knots;
}

std::vector<double> CAGD::GeomBSplineSurface::IsoVKnots() const
{
    int s = VDegree(), e = m_vknots.size() - VDegree();
    std::vector<double> knots;
    for (int i = s; i < e; i++)
        knots.push_back(m_vknots[i]);
    return knots;
}

int CAGD::GeomBSplineSurface::NbIsoU() const
{
    return m_uknots.size() - UDegree() * 2;
}

int CAGD::GeomBSplineSurface::NbIsoV() const
{
    return m_vknots.size() - VDegree() * 2;
}

std::vector<CAGD::BPnt *> CAGD::GeomBSplineSurface::MakeReference()
{
    std::vector<BPnt *> reference;
    for (int i = 0; i < m_cpts.size(); i++)
        for (int j = 0; j < m_cpts.front().size(); j++)
            reference.push_back(&m_cpts[i][j]);
    return reference;
}